#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define NUM_NEURONS	5
#define NUM_SYNAPSES	4
#define N		2
#define LEARNING_FACTOR 0.01
#define ACCEPTABLE_ERROR 0.01
#define TEMP_FACTOR	0.999
#define	NUM_SAMPLES_FOR_AVERAGE	1000.0
#define NUM_ITERATIONS_FOR_CONVERGENCE 10000


double n[NUM_NEURONS];
double b[NUM_NEURONS];
double w[NUM_SYNAPSES];

double positive_synapses_statistics[NUM_SYNAPSES];
double negative_synapses_statistics[NUM_SYNAPSES];
double positive_bias_statistics[NUM_NEURONS];
double negative_bias_statistics[NUM_NEURONS];

double data[N][2];


void
read_data()
{
	data[0][0] = 0.0; data[0][1] =  1.0;
	data[1][0] = 0.0; data[1][1] =  1.0;
}


void
print_neurons()
{
	int i;
	
	for (i = 0; i < NUM_NEURONS; i++)
		printf("%0.2lf, ", n[i]);
		
	printf("\n");
}


void
print_biases()
{
	int i;
	
	for (i = 0; i < NUM_NEURONS; i++)
		printf("%0.2lf, ", b[i]);
		
	printf("\n");
}


void
print_synapses()
{
	int i;
	
	for (i = 0; i < NUM_SYNAPSES; i++)
		printf("%0.2lf, ", w[i]);
		
	printf("\n");
}


void
randomize_neurons()
{
	int i;
	
	for (i = 0; i < NUM_NEURONS; i++)
		n[i] = (double) (rand() % 2);
}


void
randomize_bias()
{
	int i;
	
	for (i = 0; i < NUM_NEURONS; i++)
		b[i] = ((double) ((double) rand() / (double) RAND_MAX) * 2.0 - 1.0);
}


void
randomize_synapses()
{
	int i;

	for (i = 0; i < NUM_SYNAPSES; i++)
		w[i] = ((double) ((double) rand() / (double) RAND_MAX) * 2.0 - 1.0);
}


void
build_network()
{
	randomize_neurons();
	randomize_bias();
	randomize_synapses();
}


double
logit(double t, double T)
{
	return (1.0 / (1.0 + exp(-t / T)));
}


double
logit_clamped_function(int i, double T)
{
	double value;
	
	value = n[i - 1] * w[i - 1];
	value += n[i + 1] * w[i];
	value += b[i];
	
	return (logit(value, T));
}


double
logit_sample(double logit)
{
	double sample;
	
	sample = (double) rand() / (double) RAND_MAX;
	if (sample > logit)
		return 0.0;
	else
		return 1.0;
}


void
run_hidden_net(double T)
{
	int i;
	double logit;
	
	for (i = 1; i < NUM_NEURONS - 1; i++)
	{
		logit = logit_clamped_function(i, T);
		n[i] = logit_sample(logit);
	}
}


void
run_full_net(double T)
{
	int i;
	double value;
	
	run_hidden_net(T);
	
	value = n[1] * w[0] + b[0];
	n[0] = logit_sample(logit(value, T));

	value = n[NUM_NEURONS - 2] * w[NUM_SYNAPSES - 1] + b[NUM_NEURONS - 1];
	n[NUM_NEURONS - 1] = logit_sample(logit(value, T));
}


void
update_synapses_statistics(double *statistics)
{
	int i;
	
	for (i = 0; i < NUM_SYNAPSES; i++)
		statistics[i] += n[i] * n[i + 1];
}


void
update_neurons_statistics(double *statistics)
{
	int i;
	
	for (i = 0; i < NUM_NEURONS; i++)
		statistics[i] += n[i];
}


void
positive_phase()
{
	double T;
	double logit;
	int i, input_sample;
	
	for (i = 0; i < NUM_SYNAPSES; i++)
		positive_synapses_statistics[i] = 0.0;
	for (i = 0; i < NUM_NEURONS; i++)
		positive_bias_statistics[i] = 0.0;
	
	for (input_sample = 0; input_sample < N; input_sample++)
	{
		n[0] = data[input_sample][0];
		n[NUM_NEURONS - 1] = data[input_sample][1];
		
		for (T = 100000.0; T > 1.0; T = T * TEMP_FACTOR)
			run_hidden_net(T);

		for (i = 0; i < NUM_SAMPLES_FOR_AVERAGE; i++)
		{
			run_hidden_net(1.0);
			update_synapses_statistics(positive_synapses_statistics);
			update_neurons_statistics(positive_bias_statistics);
		}
	}
	for (i = 0; i < NUM_SYNAPSES; i++)
		positive_synapses_statistics[i] /= (NUM_SAMPLES_FOR_AVERAGE * (double) N);
	for (i = 0; i < NUM_NEURONS; i++)
		positive_bias_statistics[i] /= (NUM_SAMPLES_FOR_AVERAGE * (double) N);
}


void
negative_phase()
{
	double T;
	double logit;
	int i, input_sample;
	
	for (i = 0; i < NUM_SYNAPSES; i++)
		negative_synapses_statistics[i] = 0.0;
	for (i = 0; i < NUM_NEURONS; i++)
		negative_bias_statistics[i] = 0.0;
	
	for (input_sample = 0; input_sample < N; input_sample++)
	{
		for (T = 100000.0; T > 1.0; T = T * TEMP_FACTOR)
			run_full_net(T);
		
		for (i = 0; i < NUM_SAMPLES_FOR_AVERAGE; i++)
		{
			run_full_net(1.0);
			update_synapses_statistics(negative_synapses_statistics);
			update_neurons_statistics(negative_bias_statistics);
		}
	}
	for (i = 0; i < NUM_SYNAPSES; i++)
		negative_synapses_statistics[i] /= (NUM_SAMPLES_FOR_AVERAGE * (double) N);
	for (i = 0; i < NUM_NEURONS; i++)
		negative_bias_statistics[i] /= (NUM_SAMPLES_FOR_AVERAGE * (double) N);
}


void
synapse_update()
{
	int i;
	
	for (i = 0; i < NUM_SYNAPSES; i++)
		w[i] = w[i] + LEARNING_FACTOR * (positive_synapses_statistics[i] - negative_synapses_statistics[i]);
}


void
bias_update()
{
	int i;
	
	for (i = 0; i < NUM_NEURONS; i++)
		b[i] = b[i] + LEARNING_FACTOR * (positive_bias_statistics[i] - negative_bias_statistics[i]);
}


int
convergence_test()
{
	int i;
	double error = 0.0;
	
	for (i = 0; i < NUM_SYNAPSES; i++)
		error += (positive_synapses_statistics[i] - negative_synapses_statistics[i]) * (positive_synapses_statistics[i] - negative_synapses_statistics[i]);
		
	error = sqrt(error);
	fprintf(stderr, "%lf\n", error);
	
	if (error < ACCEPTABLE_ERROR)
		return (1);
	else
		return (0);
}


void
train()
{
	int converged = 0;
	int num_iteractions = 0;
	
	while (!converged && (num_iteractions < NUM_ITERATIONS_FOR_CONVERGENCE))
	{
		positive_phase();
		negative_phase();
		synapse_update();
		bias_update();
		converged = convergence_test();
		num_iteractions++;
	}
	printf("converged = %d\n\n", converged);
}


void
run()
{
	int i;
	double T;
	
	randomize_neurons();
	for (T = 100000.0; T > 1.0; T = T * TEMP_FACTOR)
		run_full_net(T);
	
	for (i = 0; i < 20; i++)
	{
		run_full_net(1.0);
		print_neurons();			
			
	}
}


int
main()
{
	srand(time(NULL));
	
	build_network();
	read_data();
	train();
	run();
	printf("\n");
	print_synapses();
	printf("\n");
	print_biases();
}
